
## "Measuring Political Support and Issue Ownership Using Endorsement Experiments,
## with Application to Militant Groups in Pakistan"
rm(list=ls())
library(R2jags)
library(foreign)

## which question to include?
q.include <- 1:4

afghan <- read.csv("../data/FullSurveyClean.csv", na.strings = ".")

## endorser
afghan$endorser <- (!is.na(afghan$Q5.1) | !is.na(afghan$Q5.2) | !is.na(afghan$Q5.3)  | !is.na(afghan$Q5.4) | !is.na(afghan$Q5.5) | !is.na(afghan$Q5.6))*1 + (!is.na(afghan$Q5.1A) | !is.na(afghan$Q5.2A) | !is.na(afghan$Q5.3A)  | !is.na(afghan$Q5.4A) | !is.na(afghan$Q5.5A) | !is.na(afghan$Q5.6A))*2 + (!is.na(afghan$Q5.1B) | !is.na(afghan$Q5.2B) | !is.na(afghan$Q5.3B)  | !is.na(afghan$Q5.4B) | !is.na(afghan$Q5.5B) | !is.na(afghan$Q5.6B))*3

## recoding Don't Know as Missing
afghan$Q5.1[afghan$Q5.1 > 5] <- NA
afghan$Q5.1A[afghan$Q5.1A > 5] <- NA
afghan$Q5.1B[afghan$Q5.1B > 5] <- NA
afghan$Q5.2[afghan$Q5.2 > 5] <- NA
afghan$Q5.2A[afghan$Q5.2A > 5] <- NA
afghan$Q5.2B[afghan$Q5.2B > 5] <- NA
afghan$Q5.3[afghan$Q5.3 > 5] <- NA
afghan$Q5.3A[afghan$Q5.3A > 5] <- NA
afghan$Q5.3B[afghan$Q5.3B > 5] <- NA
afghan$Q5.4[afghan$Q5.4 > 5] <- NA
afghan$Q5.4A[afghan$Q5.4A > 5] <- NA
afghan$Q5.4B[afghan$Q5.4B > 5] <- NA
afghan$Q5.5[afghan$Q5.5 > 5] <- NA
afghan$Q5.5A[afghan$Q5.5A > 5] <- NA
afghan$Q5.5B[afghan$Q5.5B > 5] <- NA
afghan$Q5.6[afghan$Q5.6 > 5] <- NA
afghan$Q5.6A[afghan$Q5.6A > 5] <- NA
afghan$Q5.6B[afghan$Q5.6B > 5] <- NA
  
Y <- array(NA, c(nrow(afghan),6,3))
Y[,1,1] <- afghan$Q5.1
Y[,1,2] <- afghan$Q5.1A
Y[,1,3] <- afghan$Q5.1B
Y[,2,1] <- afghan$Q5.2
Y[,2,2] <- afghan$Q5.2A
Y[,2,3] <- afghan$Q5.2B
Y[,3,1] <- afghan$Q5.3
Y[,3,2] <- afghan$Q5.3A
Y[,3,3] <- afghan$Q5.3B
Y[,4,1] <- afghan$Q5.4
Y[,4,2] <- afghan$Q5.4A
Y[,4,3] <- afghan$Q5.4B
Y[,5,1] <- afghan$Q5.5
Y[,5,2] <- afghan$Q5.5A
Y[,5,3] <- afghan$Q5.5B
Y[,6,1] <- afghan$Q5.6
Y[,6,2] <- afghan$Q5.6A
Y[,6,3] <- afghan$Q5.6B

N <- nrow(Y)
J <- length(q.include) ## number of issue areas / policies
K <- 3 ## number of treatment groups, including the control group

Y <- Y[,q.include,]

Y2 <- matrix(NA,N,J)
endorser <- matrix(afghan$endorser, N, J)
for (i in 1:N){
  for (j in 1:J){
    Y2[i,j] <- Y[i,j,endorser[i,j]]
  }
}
Y <- Y2

## drop those who are missing in all responses
#all.missing <- (apply(is.na(Y), 1, sum) == 4)
#afghan <- afghan[!all.missing,]
#Y <- Y[!all.missing,]
#endorser <- endorser[!all.missing,]
#N <- nrow(Y)

## other variables
afghan$province.id <- as.numeric(as.factor(as.character(afghan$M5)))
afghan$district.id <- as.numeric(as.factor(as.character(afghan$DISID)))
afghan$village.id <- as.numeric(as.factor(as.character(afghan$PPLID)))

village <- afghan$village.id
district <- unique(afghan[,c("village.id", "district.id")])[order(unique(afghan[,c("village.id", "district.id")])[,1]),2]
province <- unique(afghan[,c("district.id", "province.id")])[order(unique(afghan[,c("district.id", "province.id")])[,1]),2]
n.village <- length(unique(village))
n.district <- length(unique(district))
n.province <- length(unique(province))

##########################
##  The Village Model  ##
##########################
data <- list ("N", "J", "K", "Y", "endorser", 
	"village", "district", "province", "n.village", "n.district", "n.province")
inits <- function() {
  list(alpha0=matrix(seq(-2,2,length.out=(4*J)),J,4), 
       beta0=runif(J), x=rnorm(N), s=matrix(0, N, J), 
       delta.village0=rnorm(n.village, sd = 5), delta.district0 = rnorm(n.district, sd = 5),
       delta.province = rnorm(n.province, sd = 5),
       lambda.village0=matrix(c(rep(0,n.village), rep(0,n.village*(K-1))), n.village, K),
       theta.district0=matrix(c(rep(0,n.district), rnorm(n.district*(K-1), sd = 5)), n.district, K),
       theta.province0=matrix(c(rep(0,n.province), rnorm(n.province*(K-1), sd = 5)), n.province, K),
       omega=c(0,runif(K-1)), sigma=runif(n.district), psi.district=matrix(c(rep(0,n.district), runif(n.district*(K-1))), n.district, K),
       psi.province=matrix(c(rep(0,n.province), runif(n.province*(K-1))), n.province, K)
       )
}
params <- c("sd.x", "beta", "alpha", "psi.district", "psi.province", "omega",
            "sigma.district", "sigma.province", "delta.village", "delta.district",
            "delta.province", "lambda.village", "theta.district", "theta.province",
            "s", "mean.correct") 
model <- "ModelVillageDistrictProvince.R" 
fit <- jags(data, inits, params, model.file=model, n.iter=5000, n.chains = 3)

sink(paste("results/endorseResultsNoCovariates.Rout", sep =""))
print(fit)
sink()

save(fit, file = "results/endorseResultsNoCovariates.RData")

check <- any(fit$BUGSoutput$summary[, "Rhat"] > 1.2) ##Set standard for convergence
counter <- 0
while(check & counter < 25){ ##Continue until convergence is reached or 25 attempts 
	fit <- update(fit, n.iter=5000, n.thin=10)
	check <- any(fit$BUGSoutput$summary[, "Rhat"] > 1.2) ##True if any parameters have not reached convergence
	counter <- counter + 1
        sink("results/endorseResultsNoCovariates.Rout")
        print(fit)
        sink()
	save(fit, file = "results/endorseResultsNoCovariates.RData")
}  

  

